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The objective of this investigation is to study the damage mechanics of composite 


structures using a micromechanical approach for determining strength and stiffness 
degradation of the composite structures as damages, such as matrix cracking and fiber 
breakage, progress. The micromechanical cell method provides for analysis of stress at 
the fiber and matrix level while providing smeared composite properties for global 
structural analysis. As a result, the damage and failure criteria are expressed in terms of 
the fiber and matrix stress level of the composite structure. A correlation for stiffness 
reduction due to transverse cracking of a ceramic matrix composite under tensile loading 
is implemented in a three-dimensional finite element model. Next, thermal residual 
stresses from fabrication of the ceramic matrix composite are incorporated into the 
analysis. Finally, the finite element method is applied to a polymer matrix composite 
laminate with a center hole in order to study the progression of damage and final failure 
during tensile loading. The comparisons between the present predictions and the 


experimental results for the previous examples are very good. 
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L INTRODUCTION 


A. BACKGROUND 

Modem structural designs increasingly incorporate man-made composite materials 
in applications that require components with special material properties which are 
unavailable from conventional metals or alloys. From the structural mechanics viewpoint, 
composites are typically used to improve either the stiffness-to-weight ratio or the 
strength-to-weight ratio of a structural member. Naval shipboard applications of 
composite materials are typically dependent on these weight savings when a composite 
material is used to fulfill a specified stiffness or strength requirement. Current naval 
composite applications include main propulsion shafting, decks, superstructures, small 
equipment foundations and sraphitesieintorcedcesly (GRP) auxiliary system piping. 
Although the composite material itself may actually cost more than its conventional 
metallic predecessor, incorporation of lightweight composite materials usually reduces 
total ship cost due to synergistic savings in the entire ship structure because of the 
reduced total ship weight. 

Composite properties such as strength and stiffness are dependent upon the volume 
fraction of the fiber and the individual properties of the constituent fiber and matrix 
materials. The variation of the fiber volume fraction for a given composite allows the 
designer greater flexibility when incorporating composite materials irto a structure. With 


this flexibility comes added complexity due to the fact that the optimum volume fraction 





for a given stiffness may not be optimum for composite strength. And, composite 
material properties may vary due to damage accumulation such as matrix transverse 
cracking during the loading history of the composite member. This current research 
investigates the analysis of composite material properties, such as stiffness and strength, 


and their degradation using a three dimensional micromechanical model. 


B. DETERMINING COMPOSITE MATERIAL PROPERTIES 

The study of composite materials and structures can be undertaken from two 
different approaches: micromechanical and macromechanical. In the micromechanical 
approach, the properties of the constituent fiber and matrix ‘aeienalzand their interaction 
through stress-strain constitutive relations are analyzed in order to predict the overall 
behavior of the composite structural member. In the macromechanical approach, the 
properties of the constituent fiber and matrix materials are averaged or smeared to 
produce a set of pseudo-homogeneous properties for the composite structural member. 
The macromechanical approach has the advantage of requiring less detailed modeling in 
that individual fiber and matrix properties are only used initially in determining the 
smeared composite properties. Finite element models using the macromechanical 
approach are somewhat less complex and typically require less computational time than 
those models vhich are based on the micromechanical approach. However, the 
micromechanical approach provides more useful information such as the fiber and matrix 
stress and strain which are typically used for failure criteria or strength degradation 


computations. Thus, the more complex micromechanical approach is more useful when 





analyzing damage mechanics and failure of composite material structures. 


Under normal, non-damaging or elastic loading conditions, the stiffness properties 
of aligned, continuous, fiber composites can be predicted by the simple rule of mixtures 
(ROM). The ROM model uses the strength of materials approach and is based upon 
assumptions that the composite is loaded at low strain levels in the elastic region where 
no damage occurs in the fibers or matrix, that the fibers have uniform properties and are 
aligned parallel throughout the composite, that the matrix and fibers have a no-slip 
perfectly bonded interface, and that, for the longitudinally applied load case, there is equal 
strain in the fiber and matrix [Ref.1]. For composites which approximate these 
assumptions and which are loaded in the longitudinal direction along the fiber, ROM 
provides reasonable values for composite stiffness. However, for transverse loading with 
respect to the fiber axis, ROM does not accurately predict the composite stiffness. 

Micromechanical models based on the theory of elasticity provide better overall 
results for aligned, parallel fiber composites [Ref. 2]. Noteworthy among the 
micromechanical models are the relations developed by C.C. Chamis [Ref. 3]: Chamis 
assumed that the matrix is isotropic and that the fibers are orthotropic. A composite with 
such constituents would have one plane of symmetry and therefore would be transversely 
isotropic. With the condition of transverse isotropy, Chamis reduced the problem from one 
with three dimensions to one with two dimensions. This method resulted in five 
independent properties in order to define a stress-strain constitutive relation for the 
smeared composite. Chamis and Sendeckyj [Ref. 4] provide an excellent survey of the 


early micromechanics approach to predicting composite stiffness. Jones [Ref. 5] 





succinctly presents methods for bounding stiffness properties of composites using 
variational energy approaches with classical elasticity theory and also describes the 
contiguity approach used to develop the famous Halpin-Tsai relations. These previous 
methods provide generally reasonable values for composite stiffness. But, these early 
micromechanical methods are often not sufficiently accurate or readily applied to 
structural computations involving composites. And, these methods are based on two- 
dimensional, linear elastic analysis while the real composite member is three dimensional, 
and typically behaves in nonlinear manner over its useful loading range [Ref. 6 ]. 

In recent years, there has been significant study of three-dimensional 
micromechanical models. Dvorak and Bahei-E!-Din [{Ref. 7] formulated an axisymmetric 
model for an elastic-plastic constitutive relation based on cylindrical fibers with 
surrounding matrix material. The axisymmetric assumption was considered physically 
valid and reduced the complexity of the three dimensional elasticity analysis for the 
micromechanical modelling of elastic-plastic deformation. Aboudi [Ref. 8] introduced the 
micromechanical method of cells based on the assumption that the composite material 
consisted of two repetitive phases of fiber and matrix materials. By assuming the 
repetitive or periodic nature of the fiber and matrix array, Aboudi was able to simplify 
his model by using a representative cell which consisted of four subcells. Although a 
view of his model appears to desc ibe a square cell with one subcell of fiber and three 
subcells of matrix, the geometry of the fiber and matrix is not restricted due to calculation 
of the interface conditions on an average bases. Aboudi introduced the method of cells 


starting with simple unidirectional fiber composites and then extended his method to 


discontinuous short fiber composites. Aboudi's micromechanical model is mathematically 
complex and computational intensive. Others, such as Pecknold [Ref. 9] noted that 
Aboudi's model forms the basis for a finite element model. Pecknold conducted an 
investigation of a simplified unit cell model. Kwon used a micromechanics model [Ref. 
10, 11, 12] by focusing on the fiber and matrix stresses at the micromechanical level. 
Kwon later refined his original cell method [Ref. 6]. Kwon's micromechanics models are 
especially applicable for the investigation of composite damage because it considers both 
the fiber and matrix stress at the micromechanical level and thus allows specific fiber and 
matrix yield and/or failure criteria to be applied. The refined Kwon model [Ref. 6] 


formed the basis for this current work. 


C. THE MICROMECHANICAL CELL MODEL 

Similar to the Aboudi cell model, the Kwon model considers the composite as a unit 
cell composed of four subcells: one fiber subcell and three matrix subcells. The unit cell 
is represented by a three dimensional solid, the rectangular parallelpiped. Based on 


symmetry, only one-quarter of the total cell may be modeled as shown in Figure 1.1. 











Figure 1.1 - The Kwon Micromechanical Model Unit Cell And Subcells 








Note that the size of each cell is dependent on the fiber volume fraction. Kwon 
then expresses the composite stresses and strains as a function of the subcell stresses, the 


subcell strains, and the volume fraction as shown in Equations (1.1) and (1.2): 


3, = V9; +, (i -WW, kot +, ( -WW, bi o( -Wv, Jo! i,j =1,3 (1.1) 
F =Vel+/V, ( -v, ks +, (i -WW, ks + (i -\W, Je! ij =1,3 (1.2) 


where G, and &, are composite stresses and strains, o%, and e%, are subcell (a = a, b, c, 
or d) stresses and strains, and V; is the fiber volume fraction. 

Note that subcell 'a’ represents the fiber and subcells 'b’, 'c', and 'd' represent the 
matrix. Thus, the composite smeared stress or strain values are equal to the volume 
average of the subcell stresses and strains. The condition of stress continuity is satisfied 


at the subcell interfaces as expressed by Equation (1.3). 


ao_ ob ¢ od a oe b Od 
O22 =O27 , G27 =O. , G33 =G33 , 033 = 933 (1.3) 


b d b d 
Giz =Oi2 , yp =n , Fy =Gi3 , Oy =G43 

a _ ob a _oe a od 

G23 = 9x3, G23 =G33 , G23 = G25 
Each subcell may have different strains, but the following strain compatibilities are 
assumed. The longitudinal strains of the subcells are equal and the sums of the transverse 
strains of the subcell in either transverse direction are also equal. The shearing strains 
must also satisfy similar conditions. Expressed mathematically, the strain conditions are 


shown in Equation (1.4). 








a b a c a d 
Gy, =O, , Sy FS . Sn FS 
a b c d a c b d 
E22. F822 = E22 4822, F353 FE33 = E33 TE35 (1.4) 
s b _ cc a a ¢ _ ob d 
Ei2 FEy2 FSi TE.» &y3 FE 3 = Eis FEyy 
The constitutive relation for each subcell is expressed by the generalized Hooke's 


law of Equation (1.5). 


of =E%.e4 i,j,k, =1,2,3 anda =a,b,c,d (1.5) 


Solving Equations (1.1) to (1.5) together yields the expression for smeared 
composite material properties in terms of fiber and matrix material properties and their 
volume fractions. In addition, the equations provide the following sequential relations: 
global displacements —»> composite strains —> fiber and matrix strains —> fiber and 
matrix stresses —> composite stresses. 

The global displacements are obtained from the finite element analyses of fibrous 
composite structures. For damage progression studies, the calculation of stress and strain 
at the fiber and matrix level allows failure criteria to be applied at that micro-level. Thus, 
numerical modeling using this method can provide insight into the micro-level failure 
process for composite materials. 

For this current work, the Kwon micromechanical method of cells was implemented 
within the finite element program in order to perform several function:: First, the new 
micromechanical model was used to determine the smeared composite properties for use 
in computing the basic finite element stiffness matrix. Secondly, the model was 


implemented in the finite element post-processing routines to convert the finite element 














displacement results into stresses and strains at the fiber and matrix micro-level. The 
fiber and matrix stress levels were used for either stiffness reduction or the application 


of failure criteria. 











IL FINITE ELEMENT MODEL DEVELOPMENT 


Finite element numerical analysis was used to solve the three dimensional elasticity 
problem in order to determine stress and strain of the fiber, matrix or the overall 
composite specimen. The finite element method transforms the partial differential 
equations to a system of algebraic equations. The primary solution was obtained in terms 
of displacement. A solution post-processor subroutine converted the displacement solution 
into composite strain. The Kwon micromechanical method of cells provided the algorithm 
to solve for local cell strains, local cell stresses and composite stress. Stress values 
determined at the micromechanical cell !evel provided the dependent variable for stiffness 
reduction or failure relations. Details of the finite element method derivation are given 


below. 


A. ANALYTICAL DERIVATION 
The derivation for the force equilibrium equations of a three dimensional solid body 
which is experiencing negligible body forces is found in any advanced solid mechanics 


or elasticity textbook [Ref. 13, 14, 15, 16]. These equations are shown in Equation (2.1). 

















oo ,, ot 
X+_% +_™% =0 
Ox Oy a 
Ry {Py Mn og (2.1) 
Ox Oy a 
ot. oO, 0s, 
at, ee =0 
Ox Oy a 











These three equilibrium equations are written in terms on nine stress variables. However, 
only six of the stress variables are independent due to the requirement for moment 
equilibrium. Application of moment equilibrium conditions on the unit solid element 


yields Equation (2.2). 


Ty = Sx 
t, =t, (2.2) 
Vy “ty 


The unit solid element is shown in Figure 2.1 below with stress terms as indicated to 


establish a reference for the sign convention and notation. 


y 
Oy 
ty 
Tex = 


Figure 2.1 Unit Solid Cell With Stress Components 








Let ‘u', 'v', and 'w' represent displacement in the ‘x’, ‘y', and ‘z’ directions 


respectively. Then, the relationship of strain to displacement, assuming small 


displacements, can be written as shown in Equation (2.3). 


2.3 
Zr oe 
% lay ox)’ * lay a)’ * la ox 


The strain-displacement relations introduce six strain terms as unknown variables. Also, 


at this point, 'u', ‘v', 


and 'w' are unknown. These equations can be written more 


conveniently in matrix form as in Equation (2.4). 


9 9 9 
Ox 
g, 0 2 0 
: oy 
: 0 Oo ie: u 
Bee coal ee (2.4) 
0 0 
er ay 0 w 
ee 
d Oo 
Y 0° 2 
XZ Oz Oy 
Go & 
Oz Ox 


Next, consider the constitutive relationship between stress and strain of Equation (2.5). 


(0) -[D]{e] Qs) 





In Equation (2.5), {ao} is a6 x I’ column vector of stress, {e} is a ‘6 x I' column vector 
of strain, and [D] is the '6 x 6' matrix of material properties. Equations 2.1, 2.3, and 2.5 
combine to give 15 equations with 15 unknown variables. Unfortunately, the equations 
include several partial differential equations. 

The derivation of the finite element equations from the three dimensional elasticity 
equation above is based on course notes and commonly available current textbooks [Ref. 
17, 18, 19, 20]. The method of weighted residuals as modified by Galerkin is employed 
to develop the displacement based finite element equations. Eight noded isoparametric 
rectangular parallelpiped elements are used for the formulation. Linear shape functions 
are used for analytical and computational simplicity. Detailed derivation according to 
these precepts is continued below. 

The first step in converting the partial differential equations into algebraic equations 
is to apply the method of weighted residuals to the three dimensional stress equilibrium 
Equation (2.1). To this end, each of the three equations of Equation (2.1) are multiplied 
by a weight function, W; (i=1, 2, 3), which is continuous over the physical domain of the 
problem. Then, the three new product equations are integrated over the entire problem 
domain. The goal is to chose weight functions W, which are orthogonal to the initial 
residuals of the equilibrium equations such that the integral becomes equal to zero. If 


'V' is the domain volume of the problem, the weighted residual equations are now shown 


in Equation (2.6). 








(2.6) 


At this point it should be noted that boundary conditions must be specified for Equation 


(2.1). The boundary conditions may be specified as (1) essential or geometric boundary 


conditions where some surface displacements are specified, (2) natural or stress boundary 


conditions where surface tractions are specified, or (3) a combination of these types of 


boundary conditions. Before applying boundary conditions, further manipulation of the 


weighted residuals is required as shown by Equation (2.7): 


Integration by parts, 


commonly referred to as Gauss’ theorem or the divergence theorem when applied in three 


dimensions, will be performed in order to weaken the continuity requirements on the 


approximate solution. 


ee ae ee 


h= fff He + 0 tte ley + Sf (arate snes 2-7) 


1 Sf {ee Bho Ef «fF (ans SAMs #0 HAS 


In Equation (2.7), 'S' is the domain surface boundary, 'V' is the domain volume, and 'n,’, 


, and 'n,' are outward unit surface normal directions cosines in the 'x 


14 


, and 'z' 








directions, respectively. Now, define the boundary stress conditions by surface tractions 


terms as shown in Equation (2.8): 


$, =On, +t n, +Tn, 
, =n, +o.n, +t, 0, (2.8) 


o, =t “ +1, +on, 


Equation (2.7) can be written in matrix form with Equation (2.8) incorporated: 


ga Ww, 
+T — 
~ae Y & on, 
WwW, _ mM, mW, _ (2.9 
Lf} } % a — + Ey xe dV = J Pst, dS ) 
mW,  W, mM, On, 


Equation (2.9) can be further modified by separating the column vector inside the volume 
integral into a product of a '3 x 6’ matrix and a ‘6 x I‘ vector. This step, shown in 
equation (2.10), isolates the weighting function derivatives in the matrix and the stress 


terms in the vector. 


Oo 
W, ww wy, aw w,|| 
o,— + t,— +r, + —! 0 0 —_! Qo —! o 
a 7 Y a a yy a y 
0. 

pe ge Oa MN ak eg 8 ag ae de AO) 
Tae Ty Ke EY axe fy 
7%. hm! | 9 MH 4 Mh MIs, 
a Fy oa & gH & - 
xz 


Substitution of the strain-displacement Equation (2.3) into the constitutive relation 


given by Equation (2.5) provides another useful identity as shown by Equation (2.11). 


15 








du 
Ox 
ad 
oy x 
aw 
Oz 


o (2.11) 


[D] “4% 


[> 


I? ef 


+ 


V2? Vl Vy 
VP 


Equation (2.4) is now substituted into the left hand side column vector of Equation (2.11). 
Both Equation (2.10) and the modified Equation (2.11) are substituted into Equation (2.9) 


to produce Equation (2.12). 


a 
5099 
a 
0 20 
My og My MH 3 
ne So an 00 | la J, 
or Os ee OPH, 5 vt av = ff yr,tas 
Vv ean eae 5 lg 
0 0 Mg mM Hy & a 
a ey & 0 2 2 (2.12) 
ay 
a a 
a°@ 








Until this point, the elasticity relations have been manipulated to establish 
displacement as the variable for which to solve. Now, to discretize the problem for 
algebraic computational solution, assume that over a small domain each of the 
displacements can be represented by a polynomial. By discretizing a three dimensional 
domain, small-but-finite-sized volume elements can be formed. These finite elements 
become eight-noded. Each node can have three orthogonal displacements. Then, each 
element has a total of 24 degrees of freedom (dof). Calling ‘u;, 'v,, and 'w,’ the 


displacement at each node, displacements can be defined in terms of polynomial shape 


functions as follows in Equation (2.13). 


8 
u=})Hu; v=)oHy,; w=) Hw, (2.13) 


Also, the first partial derivatives of the shape functions exist as shown in Equation 


(2.14). 
a > OH, a y oH; Ow : OH; 
— = —l; = = —v; — —w 
a Fa ! a ax! a > a! 
8 8 8 
BSS Sig Me Eye We | 4014) 
WY i ® WY in & WY i ® 
fe Ti : On, Oo : HH, Ow : on, 
—_—= —u; —-F=)—v; =F —w 
a fai & > a! & » ax! 


The '3 x 1' displacement vector can be expressed as the product of a '3 x 24' shape 
function matrix and a '24 x 1' nodal displacement vector as shown in Appendix A. The 


product of the '6 x 3' partial differential matrix of Equation (2.12) and the '3 x 24' shape 


17 








function matrix are typically combined into one '6 x 24' matrix. This matrix is referred 
to as the 'B' matrix in this development and the full 'B' matrix is also shown in Appendix 
A. In shorthand notation shown in Equation (2.15), the 'B’ matrix can be partitioned into 


sub-matrices labelled 'B;', where i=1 to 8 for eight sub-matrices. 


BB) -[ 2) 8) 2) BB) BB | 15) 


Based on this shorthand notation, the sub-matrices are defined in Equation (2.16). 


OH. 
an 2 0 
Ox 
OH. 
U- *22 0 
oy 
0 0 as 
Oz 
[ B, | ie : where i =1 to 8 (2.16) 
oH, OH, 9 
Oy x 
jo 
Oz Oy 
oH, : dH, 
a ox 


The Galerkin method (or more exactly the Bubnov-Galerkin method) takes the 


weighting functions as equivalent to the shape functions as shown in Equation (2.17). 
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(2.17) 


ce ie -aie - 


Thus the Galerkin method only requires that the weighting function be continuous over 
small discrete intervals which correspond to the sides of the finite elements. Based on 
the Galerkin weighting functions shown above, the old weighting function matrix of 


Equation (2.12) can now be written in terms of the shape functions as per Equation (2.18). 


ow 
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Thus, Equation (2.12) can be revised as shown in Equation (2.19). 
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In Equation (2.19), {d} is the displacement vector which has grown from the original ‘3 


x I' vector of ‘u’, 'v’, and 'w' to a '24 x l' vector of terms 'u,, 'v;, and'w.. Also, the 
right hand side surface integral '3 x 1' vector has a vector for each term and is simply 
shorthand for a '24 x 1' vector of discretized surface boundary tractions. The original 
three partial differential equations of elasticity equilibrium have now been reduced to a 
matrix equation of size 24. The integrals in Equation (2.19) may be easily solved if 
simple forms for shape functions are chosen and if the modelled section geometry is 
relatively simple. 

If the modelled section geometry is complex, the finite element geometry can be 
made simple by employing a transform mapping to another reference space. In order to 
map ‘x’, 'y', and ‘z' coordinates for an irregular shaped element onto 'r’, 's' and 't' 
coordinates for a rectangular parallelpiped, the Jacobian transform matrix is required. The 
Jacobian transform matrix essentially scales each component of the original space to a 
new space. Thus, Equation (2.20) shows the Jacobian for an 'xyz' system mapped to an 


'rst' system: 


(2.20) 
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In terms of the shape functions and the nodal points, the Jacobian is expressed in 


Equation (2.21) on the following page. 
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However, since the finite element integral equation includes shape function partial 
derivatives with respect to the 'xyz' coordinate system, the inverse of the Jacobian is 
required as described. Let [ J ]'= [TI ], where [IT ] is a ‘3 x 3' matrix. The shape 


function derivatives with respect to the 'xyz' system are now expressed in Equation (2.22) 
p q 


with respect to the 'rst' system. 
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Equation (2.22) is used to calculate the [ B ] and [ B ]' matrices in the 'rst' system. 


Equation (2.21) is used to transform the volume differential. The finite element integral 


is transformed below as shown in Equation (2.23). 


fff TP 18) Mey (2) = fff lB (2118 11d | Mes (d) (2.23) 


The transformation to the 'rst' coordinate system results in simplified finite element 


integrals. The resultant transformed elements are termed isoparametric elements. The 
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term isoparametric refers to the equality between the degree of the equations for 
transforming ‘xyz’ into 'rst' coordinates and the degree of the shape functions for 
estimating displacement. The shape functions defined in terms of ‘rst' system variables 


are shown in Equation (2.24). 


N 


ww 


(2.24) 


~~ n a 


a — a 
t 


1 -r) (1 +s) (1 -t) 


In order to compute the integrals with a computer program, the Gauss quadrature 
of numerical integration is used. The volume integral is redefined as the triple summation 
from | to the number of integration points (NIP) of the integrand evaluated at Gauss 


integration points (r,, s, and t,) and multiplied by weighting factors as shown in Equation 


(2.25). 
NIP NIP NIP 
S{fl2 PPB sl desg (4) = FY LBL PI BIS) mM, La) 


(2.25) 


The results of the numerical integration may vary over the elements in the domain of the 
model. Each of these results can be expressed as a '24 x 24' matrix which is termed the 


elemental stiffness matrix [ K, ]. 
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Recall the surface traction boundary conditions from the right-hand side of Equation 
(2.12). The integration of the directional component of the applied stress over the 
element surface area of the applied stress is equal to the external force applied to the 
element. The result of the integration is a '24 x 1' vector { F, } which is equivalent to 
the force in the 'x’, 'y', and 'z’ direction at each of the eight element nodes. Substituting 
these resultant terms into Equation (2.24) completes the finite element derivation which 


is shown in Equation (2.26). 


[K. ]{4} ={F.} (2.26) 


The finite element method involves combining these element matrix equations with 
given boundary conditions in order to form one large system of simultaneous equations 
for numerical solution. The computer programming methodology to achieve such a 


solution is described in the next section. 


B. FORTRAN COMPUTER PROGRAM DEVELOPMENT 

Finite element numerical analysis was accomplished using FORTRAN computer 
programs. The FORTRAN programs were implemented using the modular programming 
concepts and FORTRAN subroutines developed by Akin [Ref. 21]. The various 
FORTRAN programs used for this study differed from each other in only two or three 
subroutines out of approximately 55 subroutines to the main program. 

Program housekeeping data, nodal coordinates, element-node connectivity, boundary 


conditions and material properties are read into the program from one input data file via 
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several standard Akin subroutines (for example, INPUT, INVECT, INPROP, and 
APLYBC). Efficient storage of the input data, interim calculated program data, and 
calculated results is accomplished through the use of semi-dynamic storage. Two large 
data vecters are sized in the main program to store fixed point and floating point data 
respectively. Matrix or array data is converted for storage in these column vectors. 
Computation of pointers to mark data locations in the column vectors does require some 
program overhead in that up to four subroutines (for example, PT1, PT2, PT3, and PT4A) 
are called for this purpose. However, the efficient use of computer core memory through 
the semi-dynamic storage method provides a significant advantage in the number of 
elements and nodes which can be used in the model. 

Subroutine MODEL provides overall control for the rest of the program. All 
major subroutines are called from MODEL. Subroutine ASMMDS computes the element 
stiffness matrices using Gauss quadrature integration. (ASMMDS is a locally developed 
problem dependent subroutine). Two integration points are used in each coordinate 
direction at the element nodes. This results in eight independent integrations points per 
element. ASMMDS implements the subroutines (that is, MATER and COMDMX) which 
implement the Kwon micromechanical cell method to build the composite smeared 
properties. Fiber layup angle transformations are performed in subroutine STNROT. 
Finally, ASMMDS calls subroutine STORSQ to assemble the elem::ntal stiffness matrices 
into the global stiffness matrix. 

After completing ASMMDS computations, program control returns to MODEL 


where boundary conditions are applied, modifications are made to the global force column 
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vector and the global stiffness matrix (via subroutine APLYBC), and the displacement 
solution is obtained (in subroutines FACTOR and SOLVE). FACTOR uses the 
Cholesky method to factor the system square matrix into the product of a lower triangular 
matrix and its transpose. SOLVE use Cholesky-Gauss methods of forward and back 
substitutions to complete the solution. APLYBC, FACTOR and SOLVE are standard 
Akin subroutines. 

Post-processing to calculate stress, strain, stiffness reduction or damage 
accumulation is accomplished in problem dependent subroutine MSTRES. MSTRES is 
an original subroutine developed and modified as required for this study. MSTRES re- 
computes the elemental stiffness matrices, and using the displacement solution previously 
obtained, implements the Kwon cell method to compute the micro-level cell strain from 
the composite macro-level strain. The micro-level cell strains are used to compute micro- 
level cell stresses and macro-level composite stresses. Stiffness reduction and damage 
criteria are applied based on the micro-level cell stresses calculated in this module. 
Program output is arranged from this module or from other standard Akin models as 
necessary. 

Program iteration to modify material properties is required for computation of 
damage accumulation. Necessary iteration control is designed into contro! module 


MODEL for damage studies which are discussed later. 
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IL STIFFNESS REDUCTION OF A CERAMIC MATRIX COMPOSITE 


A. OVERVIEW 

In a typical composite, the stiffness or Young's Modulus of the matrix is much less 
than that of the fiber. In polymer matrix composites (PMCs), the ratio of fiber stiffness 
to matrix stiffness may be on the order of 100. In ceramic matrix composites (CMCs), 
the ratio of fiber stiffness to matrix stiffness is usually on the order of 10 or lower. The 
ratio of fiber failure strength to matrix failure strength follows an analogous relationship 
as the stiffness rations for PMCs and CMCs. Ceramic matrix materials are typically 
stiffer, stronger and more brittle than polymer matrix materials. However, PMCs _ tend 
to fail based on the stress and strain at the fiber level, while CMCs tend to fail based on 
relatively lower stresses and strains at the siaarik level (Ref. 22]. Matrix cracking (or, 
more precisely, matrix micro-cracking) is one of the causes of initial failure of all 
composites, and, in particular, for CMCs. When the matrix develops cracks, the adjacent 
fibers must carry additional load. Thus, the effect of matrix cracking is seen on the 
macro-level as a reduction in stiffness for the composite structure. 

Much previous research had been directed to model the matrix cracking 
phenomenon. Aveston and Kelly (Ref. 23] extended the models of the late 1960': for 
understanding the stress-strain behavior of composites with matrix cracks. They used 
rigorous shear-lag analysis to calculate the fiber stresses when the matrix cracks in a 


unidirectional fiber composite plate. They implement a "strain energy versus crack 
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growth" approach to relate the interfacial fiber-matrix shear stress to the unbonded matrix 
longitudinal stress after matrix cracking has occurred. Hahn and Tsai [Ref. 25] develop 
a bilinear model to explain stress-strain behavior before and after matrix cracking. Hahn 
and Tsai developed their model with experimental results for a [0°/90°), glass fiber, epoxy 
matrix composite loaded in a uni-stress condition. Garrett and Bailey [Ref. 25] and 
Parvizi et al. [Ref. 26] tested cross-ply glass fiber, polyester matrix composite materials 
to observe matrix cracking. They validated some of the earlier theoretical work of 
Aveston and Kelly. They also studied the effect of applied stress and matrix crack 
spacing. Reifsnider and his team [Ref. 27] identified a matrix crack pattern in terms of 
a Characteristic Damage State related to applied stress and based on their experimental 
work with composite laminates with multiple fiber layup angles. Talreja [Ref. 28, 29] 
uses a complex, rigorous continuum mechanics approach to derive constitutive relations 
which essentially extends the work of the authors cited above. However, Talreja’s 
approach is based on factors adjusted by a curve-fitting to experimental results. Talreja's 
method provides some insight into the stiffness reduction from matrix cracking but 
implementation of his model requires significant experimental data. 

This research attempts to extend the prior studies by developing a correlation for 
composite stiffness reduction due to matrix cracking. The correlation is based on results 
from simple tensile testing of a composite specimen. The tensile testing provides 
characterization of stiffness reduction versus applied load. While this correlation approach 


is not predictive from first principles and basic material properties, it provides a 
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methodology to predict global structural behavior of a composite material to the micro- 


level stresses in the fiber and matrix. 


B. EXPERIMENTAL BASIS 

Wang [Ref. 30] performed comprehensive testing on continuous Nicalon (SiC) fiber 
reinforced calcium aluminosilicate (CAS) composite systems. Wang focused much 
attention on standardization of specimen design and testing methodology due to lack of 
standards for testing fiber reinforced ceramic composites. Manufacturer supplied data for 
the fiber and matrix are indicated in Table 3.1 below. 


TABLE 3.1 CHARACTERISTICS OF MATRIX AND FIBERS IN THE SiC/CAS 
COMPOSITES 


Properties CAS Matrix Nicalon, SiC 
Fibers 

Elastic Modulus 

(GPa) 


0.255 Not Provided 


Wang determined the following properties for a 36% fiber volume SiC/CAS unidirectional 











composite: Elastic Modulus ranged from 124 to 131 GPa and Poisson's ratio ranged from 
0.29 to 0.30. These values agree with results of the micromechanical model. 
Wang performed uniaxial cyclic loading of unidirectional [0] , [0,/90/0,], and 


[0,/90,/0,] composites. Stiffness reduction results of these tests are used for this study. 
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C. MATRIX STIFFNESS REDUCTION CORRELATION 

Kwon and Berner [Ref. 31] proposed the following method to account for the effect 
of matrix cracking on the composite stiffness. The degradation of the matrix elastic 
modulus is assumed to be a function of the micro-level matrix stress as shown in 


Equation (3.1). 


E* =E° f (7) (3.1) 
In Equation (3.1), the superscript 'm' represents a matrix property, the subscript 'o' 
represents the original undamaged elastic modulus, and 'f is a stiffness degradation 
function which depends on the matrix cell stresses. A Weibull type distribution function 


of the matrix cell von Mises stress values is used to model the stiffness degradation of 


the matrix. The stiffness degradation Weibull function is given in Equation (3.2). 


Orme “On ifo™_ > o" 
f= exp Y oo IE Ovme ys (3.2) 
1 if Ovne S Fy 


In Equation (3.2), O",me is the equivalent von Mises stress of the matrix, o”,, is the yield 


vme 
or failure strength of the matrix, and y and B are constants of the material which are 
determined from experiment. 

A one-eighth model of a tensile test specimen was developed for input into an FEM 
program developed as described in Chapter II. The composite model consisted of two 


layers of 20 elements for a total of 40 elements. Material properties from Table 3.1, CAS 


matrix Poisson's ratio of 0.31, and fiber volume of 36% were input for the analysis. A 
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CAS matrix yield stress value of 130 MPa was used for the stiffness reduction function. 


Boundary conditions of known displacements which were derived from the Wang 
experimental strain and composite stiffness reduction data were also input into the FEM 
program. The FEM postprocessor determined the matrix and fiber cell stresses. The 
matrix von Mises level stress values were calculated at the cell level and then averaged 
throughout the model. Average matrix cell level von Mises stress values were considered 
valid as the specimen was modelled for uniaxial displacement loading. A separate 
FORTRAN program based on the Kwon micromechanical model as implemented by 
subroutine COMDMxX was run to correlate the matrix stiffness reduction to the composite 
stiffness reduction. This step essentially was directed at determining what the required 
level of matrix stiffness degradation would be for a given composite stiffness degradation. 
At this point, necessary information to perform a correlation was accumulated. This data 
is displayed in Table 3.2. Next, Equation (3.2) was manipulated (using logarithms and 
basic algebra) into the format required to perform a least squares fit for constants y and 
B. The least square fit constants were used to calculate the reduced matrix stiffness in 
order to finally compute the new predicted, reduced composite stiffness. This predicted, 
reduced composite stiffness was correlated to the experimentally determined, reduced 
composite stiffness. The least square fit constants were adjusted to optimize the 
maximum error and the root mean square error between the predicted and experimental, 
reduced composite stiffness values. These steps were performed using a Microsoft Excel 


4.0 spreadsheet, its built-in functions, and an interactive "macro" or program. 
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TABLE 3.2 STIFFNESS REDUCTION CURVE FIT DATA 


Strain % Stiffness E, Multiplier Required E,, Matrix von 
Reduction Multiplier Mises Stress 

(MPa) 
(experimental) | (experimental) (calculated) (computed) (computed) 


12.5% 0.875 


18.0% 0.820 





From the above data and analysis, y was determined to equal 0.307 and B was determined 
to equal 1.16. These material constants were assumed to be independent of the composite 
layer configuration, the material properties of the fiber, and the volume fraction of the 


fiber. 


D. FEM STIFFNESS REDUCTION RESULTS 

The correlation for stiffness reduction discussed above was incorporated into the 
computer FEM program. To achieve stiffness reduction, the program was run once with 
initial undamaged material properties to determine the matrix micromechanical stress 


values. The matrix stiffness reduction factor was calculated based on these stress values. 
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The matrix properties were degraded accordingly, and the FEM program was re-run to 


determine the final specimen displacement and strain with the matrix cracking damage 
present. Two different laminates were studied: Two cross-ply composite laminates made 
of SiC/CAS with a fiber volume of 36%. Numerical results were compared with Wang's 
experimental results as shown in Figures 3.1, 3.2, and 3.3. The numerical method for the 
cross-ply [0,/90/0,] composite predicts slightly lower stiffness reduction at low composite 
strains and slightly higher stiffness reduction at medium-to-high composite strains than 
that obtained by experiment. For the cross-ply [0,/90,/0,] composite, the numerical 
method significantly under-predicts the stiffness reduction as compared to the 
experimental results. In the latter case, post-experiment composite analysis indicated 
some fiber breakage in addition to the matrix cracking phenomenon. The fiber breakage 
offers a partial explanation for the higher than predicted stiffness reduction. Also, other 
non-modelled micro-level stresses or load transfer mechanisms might account for the 
additional stiffness reduction. 

To gain further insight into the stiffness reduction mechanism for the cross-ply 
composite laminates, the percent stiffness reduction was plotted for matrix cracking in 
each layer and for matrix cracking in the whole laminate. Figures 3.4 and 3.5 display 
these plots. One might expect a greater degree of matrix cracking to occur in the 90° 
layer, as its matrix supports a majority of the stress which results from the applied 
longitudinal load than does its fibers. But, for the [0,/90/0,] composite, most stiffness 
reduction occurred due to matrix cracks in the 0° layer. The reason is that the 0° layer 


is much thicker than the 90° layer. As a result, the 0° layer carries more load compared 
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to the 90° layer. For the [0,/90,/0,] composite, the numerical method indicates that 
stiffness reduction due to matrix cracking in both layers is comparable. The relatively 
thick 90° layer in [0,/90,/0,} yielded more contribution to the stiffness degradation, due 
to matrix cracking in the 90° layer, than the relatively thin 90° layer in [0,/90/0,]. 
Figures 3.6, 3.7, and 3.8 display the longitudinal fiber stress changes as the matrix 
cracking is modelled. For all composite laminates considered, as the matrix cracks, the 
fiber stress is increased due to load transfer from damaged matrix to the relatively stiff 
longitudinal fibers. For the cross-ply [0,/90/0,] laminate of Figure 3.7, matrix cracking 
in the 0° layer accounts for most of the stress increase in the longitudinal fibers. For the 
cross-ply [0,/90,/0,] laminate of Figure 3.8, matrix cracking in either layer contributes 
approximately the same relative share of longitudinal fiber stress increase. These fiber 


stress graphical results further indicate the utility of the Kwon cell method. 
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Figure 3.1 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A Unidirectional SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.2 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A [0,/90/0,] SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.3 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A [0,/90,/0,] SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.4 Stiffness Reduction Due To Matrix Cracking In Different Layers Of A 
[0,/90/0,} SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.5 Stiffness Reduction Due To Matrix Cracking In Different Layers Of A 
(0,/90,/0,] SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.6 Effects On Longitudinal Fiber Stress Of Stiffness Reduction Due To Matrix 
Cracking Of A Unidirectional SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 3.7 Effects On Longitudinal Fiber Stress Of Stiffness Reduction Due To Matrix 
Cracking In Different Layers Of A [0,/90/0,] SiC/CAS Composite Laminate Loaded 
Uniaxially. 
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Figure 3.8 Effects On Longitudinal Fiber Stress Of Stiffness Reduction Due To Matrix 
Cracking In Different Layers Of A [0,/90,/0,] SiC/CAS Composite Laminate Loaded 
Uniaxially. 





E. PARAMETRIC STUDY 


A parametric study for different fiber volumes in a unidirectional SiC/CAS 
composite laminate was conducted and reported by Kwon and Berner [Ref. 32]. Results 
for stiffness reduction sensitivity to the variation of fiber volume are indicated in Figure 
3.9. As the magnitude of the applied strain boundary condition was increased, the fiber 
stress also increased. This resulted in an increased matrix stress too. When matrix stress 
increased, the stiffness reduction factor also increased. Thus, composite stiffness 
reduction increased for all fiber volumes modelled when the strain was increased. At 
higher fiber volumes, the numerical model indicated smaller stiffness reduction. This 
result is expected: As the amount of fibers was increased, the share of load carried by 
the fibers increased likewise. Thus, local matrix stress and subsequent matrix cracking 
was decreased. 

Figure 3.10 indicates the longitudinal fiber stress versus strain for several different 
fiber volume fractions. The fiber stress increases with strain due to increasing load 
transfer from the damaged matrix to the undamaged fibers. As the effect of matrix 
cracking is reduced with higher fiber volumes, the rate of increase of fiber stress goes 
down. 

The parametric study results appcar reasonable. However, comparison with 
experimental data is warranted to fully validate the results. Also, experimental 
verification could substantiate the assumption that the stiffness reduction correlation, 


developed in section 'C' above, is truly independent of the fiber volume. 
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Figure 3.9 Stiffness Reduction Versus Strain For Different Fiber Volumes For A 
Unidirectional SiC/CAS Composite Laminate. 
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Figure 3.10 Longitudinal Fiber Stress Versus Strain For Different Fiber Volumes For A 
Unidirectional SiC/CAS Composite Laminate. 
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IV. STIFFNESS REDUCTION MODIFIED BY THERMAL RESIDUAL STRESS 


A. BACKGROUND 

The mismatch between the coefficient of thermal expansion (CTE) for the fiber and 
the matrix of ceramic matrix composites (CMCs) and the relatively high fabrication 
temperatures for CMCs can cause significant residual thermal stresses. These residual 
thermal stresses may have a significant effect on the performance of CMCs. 

The micromechanical matrix stress was the primary factor in the stiffness reduction 
correlation developed for the SiC/CAS as discussed in the previous chapter. The Weibull 
type correlation has a reasonable physical basis in using the ratio of the difference of 
matrix cell stress and matrix yield stress to the matrix yield stress: This ratio indicates 
some sort of threshold matrix cell stress level at which the matrix begins to develop 
cracks. Residual thermal stresses in the matrix could affect the threshold stress for the 
onset of cracking. 

A calculation was performed to estimate the residual thermal stresses in the 
SiC/CAS composite laminates which were analyzed in Chapter II. The calculated 
thermal residual stress terms were incorporated into the FEM program for the analysis of 
stiffness reduction. Results of the numerically predicted stiffness reduction were 


compared with experimental data. 
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B. DERIVATION 

Two separate calculations were performed for the thermal residual stress in the 
matrix of the SiC/CAS composites: unidirectional laminate calculations and cross-ply 
laminate calculations. The coefficient of thermal expansion (CTE) for the SiC fiber, a, 
was taken as 4.0 x 10° per degree Celsius; the CTE for the CAS matrix, a, 5.0 x 10° 
per degree Celsius. This data is based on Wang [Ref. 30]. Other basic material 
properties listed previously in Table 3.1 were also used. A temperature delta of 300°C 
was used as a typical value for the difference between fiber/matrix bond setting 
fabrication temperature to application temperature. The approximate computation 
described below relies on the Chamis relations (Ref. 3] to determine smeared properties 
for the composite laminates. A method to numerically compute the thermal residual 


stresses using the Kwon cell method is briefly outlined in Appendix B. 


1. Unidirectional Case 
Consider a two dimensional view of tl.2 fiber and matrix interface as two solid 
blocks which are adjacent to one another. A force balance must exist between the bonded 
fiber and matrix after cooldown from bond setting during fabrication as shown in 
Equation (4.1). 
o, A, =o, A (4.1) 
If each side of Equation (4.1) is multiplied by the composite thickness and then divided 


by the total composite volume, Equation (4.2) is obtained. 


o, V, =6, V (4.2) 
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The mechanical strain of the fiber and matrix is given by the basic Hooke's Law 


relationship and is equal to the respective elastic modulus divided by the respective stress 
term. The thermal strain is equal to the appropriate CTE times the temperature delta. 
Equation (4.3) shows that the sum of the mechanical strains of the fiber and the matrix 
must be equal to the magnitude of the difference of the thermal strains of the matrix and 


the fiber. 


o So, 
E ‘ =| a, a, | AT (4.3) 


When Equations (4.2) and (4.3) are combined algebraically, Equation (4.4) results. 


— | a, -a, | AT 
V (4.4) 





Although developed independently, Equation (4.4) is equivalent to the result derived by 
Sambell, et. al., [Ref. 33]. When the appropriate values are substituted into Equation 


(4.4), the resultant residual thermal stress in the matrix is approximately 15 MPa. 


2. + Cross-Ply Case 
The cross-ply development for thermal residual stress is analogous to the 
unidirectional development, but necessarily more complex due to the consideration of the 
two laminae of different fib :r layup angle and the two stress directions. In the following 
derivation, subscript 'L' represents the global longitudinal direction, subscript 'T’ represents 


the global transverse direction, superscript '0' represents the 0° fiber layup lamina, and 


42 











superscript '90' represents the 90° fiber layup lamina. Four different strain terms are 


represented in Equation (4.5). 


o o 
. = = Vir ae +a AT 
EO OEY 
90 90 
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er = Vit =p Pr AT 
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Because the strains in the global transverse and longitudinal directions must be equal, 
Equation (4.5) actually represents two independent equations with four unknowns. To add 
two required equations to the system, Equation (4.6) describes the force balance between 
the laminae. 
A’o, +A”%o; =0 (46) 
A’at +A%ar =0 
To solve for the longitudinal and transverse CTE values and elastic modulus values for 
the fiber, the equation of Chamis [Ref. 3] were used. The matrix CTE was assumed to 
be isotropic. Noting that the ratio of the areas for the 0° lamina to the 90° lamina Is two 
for the [0,/90,/0,] composite and six for the [0,/90/0,] composite, the system of four 
equations and four unknowns given by Equations (4.5) and (4.6) is solved to give residual 


thermal stress values. Matrix thermal residual stress values are shown in Table 4.1. 
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TABLE 4.1 MATRIX THERMAL RESIDUAL STRESSES 


[0,/90,/0,] Laminate [0,/90/0,) Laminate 
Matrix Thermal Residual | Matrix Thermal Residual 
Stress (MPa) Stress (MPa) 





C. FEM STIFFNESS REDUCTION WITH THERMAL RESIDUAL STRESS 

The Weibull stiffness reduction function constants were re-computed using the 
matrix residual stress values. The new y was 0.238 and the new B was 1.42. Figure 4.1 
shows the graph for the stiffness reduction of the unidirectional composite laminate. The 
new predicted stiffness reduction values show good agreement with the experimental 
results, but the overall comparison indicates mixed results for this case. The new 
predicted stiffness reduction values show slight improvement over the previous (without 
thermal residual stress) predicted values at low strains. The previous model shows 
slightly better agreement with experimental results at medium to high strain rates. Both 
predictions are within the same range of accuracy for the constants and assumptions 
employed in deriving the model. 

Figure 4.2 shows the graph of stiffness reduction for the cross-ply [0,/90/0,]} 


composite laminate. The new predicted stiffness reduction values show vastly improved 





agreement with experimental results as compared to the previous model. For medium to 
high strains, the new predicted values are extremely close to the experimental values. For 
low strain levels the predicted stiffness reduction due to matrix cracking is smaller than 
both the experimental results and the previous results. However, the shape of the stiffness 
reduction curve shows upward curvature and then levels off just short of achieving 
downward curvature. This type of behavior might occur if weak fibers in the composite 
laminate were cracked or damaged at a low strain level. From a statistical viewpoint, 
some early weak fiber breakage is a real probability. 

Figure 4.3 shows the graph of stiffness reduction with thermal residual stress for the 
cross-ply [0,/90,/0,} composite. The numerical prediction with thermal residual stress 
incorporated shows slight improvement in matching the experimental results as the 
previous numerical method. These results indicate the need for further refinement of the 
stiffness reduction correlation for thermal residual stresses and their effect on matrix 


cracking. 











Figure 4.1 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A Unidirectional SiC/CAS Composite Laminate Loaded Uniaxially. 


Figure 4.2 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A [0,/90/0,] SiC/CAS Composite Laminate Loaded Uniaxially. 
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Figure 4.3 Experimental And Predicted Stiffness Reduction Due To Matrix Cracking For 
A [0,/90,/0,] SiC/CAS Composite Laminate Loaded Uniaxially. 
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V. DAMAGE PROGRESSION AND FAILURE IN A POLYMER MATRIX 


COMPOSITE 


A. OVERVIEW 

A finite element analysis was also performed concerning the effects of flaws (in this 
case, a circular hole) on the failure behavior of composite laminates. From classical 
elasticity based on isotropic materials, an "infinite" plate containing a small circular hole 
and loaded in uniaxial tension, "o,," will experience a stress concentration factor (SCF) 
at the edge of the hole. The SCF produces a stress of magnitude 3 times 6, in a region 
outside the hole and perpendicular to the applied stress. Konisch and Whitney [Ref. 34] 
noted that the classical SCF approach Srsdueed anomalous results when applied to 
composite laminates with holes. They modified a solution to the elasticity equations for 
an isotropic material for application to a composite plate with a hole. Then, using a semi- 
inverse solution technique, they determined the stresses in the vicinity of the hole. Their 
semi-inverse solution produced a slightly improved correlation with experimental data. 
However, overall force equilibrium was not satisfied throughout the composite by this 
type of semi-inverse solution. 

The study of stress around holes in composite plates is predicated by the need to 
use bolted or pinned joints to connect such plates within structures. Bolt hole analysis 


has added complexity due to the loading in the hole area. Chang, Scott and Springer 
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{Ref. 35, 36] studied the failure of such joints in composite materials. They developed 
models for the bolt hole boundary conditions and implemented a variation of the Yamada 
and Sun [Ref. 37] failure criteria) The Yamada-Sun criteria is a quadratic equation for 
composite laminate failure prediction based on average laminate properties. More 
recently, Chang and Chang have extended the early work to study the failure of composite 
plates with non-loaded holes [Ref. 38]. The recent Chang research is based upon FEM 
analysis using composite laminate theory. Chang used statistical Weibull functions to 
determine degraded stiffness properties for the composite laminate when composite 
damage criteria are exceeded. Chang implemented a fiber damage zone failure method 
covered by Hahn and Tsai [Ref. 39]. The Chang numerical predictions show good 
agreement with experimental results. This current work attempts to similarly predict 
failure progression using the Kwon micromechanical model with its strong physical basis 
in the micromechanical properties of the composite constituents. Further finite element 


analysis and modelling remain to fully compare this present method with Chang's results. 


B. FAILURE CRITERIA AND MODELLING 

Chang implemented failure criteria at the composite level using modifications of the 
quadratic relations of Yamada and Sun. Since material failure has physical origins at the 
micromechanical level, it is desirable to establish failure criteria on material properties at 
that level. Aboudi [Ref. 8] has proposed the criteria presented in Equations (5.1) and 


(5.2). 
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| oF, | < x! (5.1) 


=| =| ai (52) 
In Equation (5.1), o°,, is the local fiber longitudinal stress and X, is the ultimate tensile 
strength of the fiber in the longitudinal direction. In Equation (5.2), o”,, is the local 
transverse stress of the matrix, o”,, is the local matrix shear stress in the longitudinal- 
transverse plane, X,, is the ultimate tensile strength of the matrix, and S,, is the ultimate 
shear strength of the matrix. Thus, equations (5.1) and (5.2) are based on local cell stress 
levels of the fiber and matrix and on the ultimate strengths of the fiber and matrix. Local 
fiber and matrix cracking damage, indicating local laminate failure, is assumed to occur 
when these criteria are exceeded. 

When the composite level stresses epic his global failure criteria, Chang 
modelled failure using a different method for both the fiber and composite. For the fiber, 
Chang degraded fiber properties by reducing the fiber transverse elastic modulus to zero 
and the fiber longitudinal elastic modulus and shear modulus according to an 
experimentally determined Weibull distribution. Chang degraded the matrix properties 
by reducing the matrix transverse elastic modulus to zero. For this work, when the fiber 
or matrix criteria of equations (5.1) and (5.2) was exceeded, the respective elastic moduli 
and shear moduli for all directional combinations were reduced by a factor to near zero. 
Actual reduction to zero was not possible using numerical methods due to resultant 


singularities in the matrix equations. 
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C. FEM DAMAGE PROGRESSION STUDY 


As noted above, comparison to and extension of the Chang research was a desired 
goal. Thus, the material chosen for this numerical analysis was based on the Fiberite 
T300/1034-C graphite epoxy composite which was also used for the Chang work. A 
cross-ply polymer matrix composite was studied: a [0/90], laminate. The fiber and matrix 
properties were "reverse engineered" from the smeared composite properties presented by 
Chang. The fiber and matrix properties were determined by an interactive trial and error 
method utilizing the Kwon cell model. Table 5.1 below shows the Chang smeared 
properties and the calculated fiber and matrix properties. Fiber volume was calculated 
to be 50%. 


TABLE 5.1 GRAPHITE EPOXY T300/1034-C PROPERTIES 


Property Composite - —- Fiber Matrix 
Chang Smeared Calculated Calculated 


Longitudinal 
Elastic Modulus 
(Msi) 


Longitudinal 
Tensile Strength 
(Ksi) 


Shear Strength 19.4 NA 23 
(Ksi) 
Transverse Tensile NA 5.0 
Strength 


















The FEM program discussed earlier in this report was modified to simulate damage 
progression. The composite laminate was simulated with a one-eighth model due to 
symmetry. The one-eighth model laminate plate was 0.5 inches wide, 4.0 inches long, 
and 0.125 inches thick. Hole radius was one-fourth of the laminate model width at 0.125 
inches. Figure 5.1 displays a scaled version of the specimen one-eighth model. A 
program control loop was implemented to increase the applied end tension up to 72.0 Ksi 
in increments of either 6.0 Ksi or 7.2 Ksi, depending on the number of load iterations 
chosen. The tension was applied in the 'x' direction. Within each load increment, a 
conditional program loop tested for failure; recorded and wrote failure output data; and, 
iterated to re-compute degraded stiffness properties, re-solve for new displacements, and 
re-calculate to determine cell stresses until the failure criteria indicated no further failures 
for the current load increment. 

For the [0/90], composite laminate, total matrix failure occurred in the 90° layer at 
approximately 36 Ksi. The matrix failure occurred perpendicular to the applied load 
exactly along the region of maximum stress according to classical theory. Total specimen 
failure occurred at 54 Ksi when the fiber in the 0° layer failed across the specimen. The 
fiber failure around the hole in the 0° layer initiated 78.5° to 90° off the axis of the 
applied load, and then progressed along at 66° to 78.5° off the axis of the applied load. 
At failure, the fiber failed simultaneously along the outer specimen elements at 66° to 
78.5° off the axis of the applied load. These results are displayed in Figures 5.2 through 
5.10 on the following pages. Chang modelled a [(0/90),], laminate for numerical 


prediction of failure. For this similar laminate with the same width-to-diameter ratio 
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(W/D=4.0), Chang observed failure at approximately 52 Ksi with the failure path about 
75° off the axis of the applied load. Thus, the current results show good agreement with 
previous research in this area. 

The [0/90}, composite laminate was re-modelled with end displacements applied as 
compared to above, when end forces which are equivalent to end stress tractions, were 
applied. Failure load was within 5% of the results above. The failure path was also 
similar to the results discussed above. 

These results differ only slightly from the Chang analytical results. A primary 
cause for the slight difference probably lies with the "reverse engineered" fiber and matrix 
properties and fiber volume used for this study. A secondary difference is the failure 
criteria applied. The failure degradation model used for this current study was somewhat 
severe, in that properties were reduced to the threshold of numerical singularity. 


However, the current failure model is desirable from the point of simplicity. 
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Figure 5.1 - X-Y Plane View Of The FEM One-Eighth Model Of A Composite Plate 
With A Circular Center Hole. 
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Figure 5.2 - Damage Progress At 12.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 


Note: The small circles indicate fiber failure. The small boxes indicate matrix failure. 
Hollow symbols indicate the 0 degree layer. Solid symbols indicate the 90 degree 
layer. 








Figure 5.4 - Damage Progress At 24.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 











Figure 5.5 - Damage Progress At 30.0 Ksi Applied Stress For {0/90}, 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 








Figure 5.6 - Damage Progress At 36.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. : 








Figure 5.7 - Damage Progress At 42.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 
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Figure 5.8 - Damage Progress At 48.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 
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Figure 5.9 - Damage Progress At 54.0 Ksi Applied Stress For [0/90], 
T300/1034-C Graphite Epoxy Composite Laminate Plate With Center 
Hole. 
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VL SUMMARY AND CONCLUSIONS 


A. SUMMARY 

The new micromechanical model for composite properties developed by Kwon [Ref. 
12] was readily incorporated into finite element analysis of composite structural 
specimens. The Kwon micromechanical model also provided for useful micro-level fiber 
and matrix stresses in the finite element post-processor subroutines. The fiber and matrix 
stresses were used to initiate stiffness reduction correlations for material property 
degradation and failure criteria for damage progression until failure. 

The stiffness reduction correlation developed for SiC/CAS ceramic matrix 
composites was implemented in the finite element analysis based on micro-level matrix 
stresses determined from the Kwon ‘aicroccnadical cell model. A stiffness reduction 
correlation was developed by using finite element analysis of a unidirectional SiC/CAS 
composite. The stiffness reduction correlation was based on a Weibull type function. 
Stiffness reduction predictions for the [0,/90/0,] SiC/CAS composite showed very good 
agreement with experimental results of Wang [Ref. 30]. Stiffness reduction predictions 
for the [0,/90,/0,] SiC/CAS showed good agreement with experimental results. 

Thermal residual stresses were incorporated into the stiffness reduction correlation 
and finite element analysis for three SiC/CAS composite laminates. The unidirectional 
laminate was used to determine the initial stiffness reduction with thermal residual stress 


correlation. The finite element analysis predictions for the [0,/90/0,} SiC/CAS composite 
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showed markedly improved agreement with experimental results. However, the revised 
prediction for stiffness reduction of a [0,/90,/0,] composite showed little discernable 
improvement in matching experimental results as compared with the stiffness reduction 
correlation without thermal residuals stress included. 

Damage progression and failure studies were performed using finite element analysis 
with the Kwon micromechanical model for a graphite/epoxy (T300/1034-C) cross-ply 
composite laminate with a center hole. The finite element analysis for this study 
predicted the failure load and path with very good correlation to the previous work of 


Chang [Ref. 38]. 


B. CONCLUSIONS 

The micromechanical cell model developed by Kwon [Ref. 12] can be readily 
implemented into finite element analysis of composite structures. The initial studies 
performed in this current work indicate good correlation with experimental work and 
previous research by other authors. Of course, further analysis of other continuous fiber 
composite laminate structures with comparison to experimental results is warranted. Also, 
the micromechanical model need to be modified for discontinuous fiber composites of 
whisker composites. 

The stiffness reduction correlation developed in Chapter III needs further r>finement. 
The predictions for the [0,/90,/0,] cross-ply laminate showed only fair correlation with 


experimental results. Incorporation of the characteristic damage state or zone concept 
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based on matrix cracking density into the stiffness reduction weibull function would 
provide a more clear relationship with the physical mechanism of the matrix cracking. 

Including thermal residual stresses in the stiffness reduction correlation generally 
improved the agreement of the predictions with experimental results. To further improve 
these predictions, a full three-dimensional thermal residual stress analysis based on the 
Kwon micromechanical model should be employed. Also, more detailed information on 
composite fabrication practice with regard to heating and temperature might be useful to 
improve the computation of thermal residual stress. 

The computation of fiber and matrix stress at the micromechanical level is an 
exceptionally powerful result of the Kwon micromechanical cell method. This method 
provides an outstanding physical foundation for future study of composite structural 


stiffness degradation, damage and failure analysis. 
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APPENDIX B - THERMAL STRESS DERIVATION USING KWON CELL 
METHOD 


The following derivation has been developed by Dr. Young W. Kwon. The 
derivation is provided here for completeness and reference. 
The foundation for the Kwon micromechanical cell method is based on the relations 


given in Equation (B.1). 
ay = Veg + VTA - (hog + WA - (Woy +(t- Woy ai = 13 
Meg Wl - eg Wl -Wegt(t- Wey W138 wn) 
let a = 5 B= Wi(t- Ys v= (1- WP 
Based on the unit cell as displayed in Figure 1.1, the compatibility equations for stress 


and strain are given in Equations (B.2) and (B.3). 


_~? c  _ od 2 _ ue b ed 
O22 =Ox2 , G22 =Ox. , F833 =933 , 033 = O53 
db ¢ d a c b d 

Bip FOy2 » Gyp =Gpp , G3 =Gjy3 , Gy3 =Ty3 (B.2) 


_ a? es os a _ od 
Oy, =923 , G23 =923 , S23 = 523 


a b ry c a ad 
Sy, =Gy > Sy =F » Sy FEy 
a b _ ve d a ¢ _ ob d 3 
E22 FE2. = E22 E22 =, 33 E33 = E33 FE33 (B.3) 


a b _Q ac d s c ob d 
Giz FE;2 = E12 FE. » yy FE,3 = G3 FE} 


The total strain is equal to the sum of the mechanical strain and the thermal strain 


as given by Equation (B.4). 
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e; = ncn), *Crermsi), (B.4) 


The generalized Hooke's law relates the mechanical stress to the mechanical strain. From 
Equation (B.4), the mechanical strain can be composed of the difference of the total strain 


and the thermal strain. Equation (B.5) shows this result. 


OF = EiuEnen) = Ene ij -( Coca), = Se Ener), (B.5) 


An inverse relationship for the mechanical stress can also be developed as shown in the 


algebraic derivation in Equation (B.6) for an orthotropic material. 
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Equation (B.7) re-writes the stress compatibility equations in terms of mechanical strain 


and thermal stress that were developed in equation (B.6). 
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E3,€), +E: +E 363; -©);, 


ee eee b 
= E28), +Ep2822 +Ep:€33 -(6, 


Pee ye ae Ok eee ae 
= E,,6,,; tE2€2, +E,€33 - 


Equation (B.8) re-writes the strain compatibility equations in a similar manner. 


a b 

1, =&) = 

a b 
2. €), =&); 


€. nN *E), =a), *e), 


3. 8, =8) 
€.J, te, 
4. €) =€) 


Ca), “Eh 


(B.7) 
=E3€1 +E5,82, +E3,833 -@), 
= E81 +E5,€ 72 +E3,6533 -(,),, 
= eh, =€, 
En), +e) 
(B.8) 


d 


Em n *e),, 


a bo oc d 
5. Ey, t8x2 =&z) Ez 


Ea), 


a c 
6. €33 +€33 


Cals *E);; 


+e), ‘ 


€n),, +e) 7 En), +e), +E), *e), 


b od 
= &33 +&33 


Ga), *E),, = Cabs Ps *Eahs *C)5, 


75 











Combining terms from Equations (B.1), (B.7), and (B.8) yields a system of equations 


where the mechanical cell strains can be expressed in terms of the composite average 


strains as shown in Equation (B.9). 


r b ¢ 4 
E22 +Ex, ~€22 —E2, =0 
a b c d 
£3; E33 +833 “Ey, = 0 
a b ¢ d _ 
ae,, +Pe,, tPe,, t¥E,, =F,, 
a b c d = 
ae,, +Bs,, +Pe,, tye, =€ 
aos bb aoa bob _ {pb a b Bo 
E2822 “E2262, tEz3€3; —E2;€33 = (E}, -EXE,, +0), -(6,) ‘Be? 
aoa coc aus cc Tfpe _pe _Ipw \C 
E3282, ~E3,€2, +E3;€3; ~E3;€35 = (E5, Eye +0), Gy; 
coc did coe dod _f{p¢d onc co d 
E822 ~E2,€2. +E2;€3; ~E23€33 = (Ei, EE, +o), 6), 


bb cdid pbb did _[pa _pd b _/~\d 
E4282. ~E32€2) +E;;€3;; ~E33€33 = (B5, EE, +(0,) (,) 


Equation (B.9) can be written in matrix/vector equation format as in Equation (B.10). 


(} Lai -(469) i 


Displacement based FEM programs will solve for composite macro-level 
displacements. Postprocessing modules first convert these displacement values into 
composite strains. Inverting Equation (B.10) a!tows solving for the micro-level cell 
strains. In summary, the solution procedure for thermal stresses becomes: 


1. Compute the smeared property matrix [D] and the thermal stress vector {9,}. 

2. Compute the displacement vector {d} from {K]{d}={F}, where vector {F} 
includes both mechanical and thermal loads. 

3. | Compute the composite strain vector{e}. 

4. Compute the fiber and matrix cell strains{e*}. 

5. Compute the fiber and matrix mechanical strains {e,"} from the equation 
{e,°} + {e,"} = {e"}. 

6. Compute the fiber and matrix cell stresses {0°} from {o*} = [E] {,*}. 

7. Compute the composite stress {6} from {o*}. 
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